Add a function to save a slice in MDH format #149

granrothge · 2023-08-25T13:09:23Z

This function will take the output of a spectra and create a slice in an MDH format.
Note the "Instrument" name is hard coded as a valid instrument is needed for Mantid to read the file. I'd prefer this to be "SpinW", but since that is not a valid Mantid instrument it can't be read by Mantid. I would like to see this requirement relaxed in Mantid, but that is a different discussion.

codecov-commenter · 2023-08-25T13:19:35Z

Codecov Report

Patch coverage: 21.89% and project coverage change: +0.22% 🎉

Comparison is base (415d143) 38.67% compared to head (32cb8b7) 38.90%.
Report is 4 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #149      +/-   ##
==========================================
+ Coverage   38.67%   38.90%   +0.22%     
==========================================
  Files         239      240       +1     
  Lines       15829    15971     +142     
==========================================
+ Hits         6122     6213      +91     
- Misses       9707     9758      +51

Files Changed	Coverage Δ
swfiles/private/sw_issymspec.m	`54.54% <0.00%> (-20.46%)`	⬇️
swfiles/sw_instrument.m	`2.36% <0.00%> (-0.16%)`	⬇️
swfiles/sw_spec2MDHisto.m	`0.00% <0.00%> (ø)`
swfiles/sw_egrid.m	`85.71% <100.00%> (+22.64%)`	⬆️

... and 2 files with indirect coverage changes

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

mducle

Code looks ok. But there are a couple of issues:

I get an error when the output file already exists - could you either document this or automatically remove / overwrite the file? (I think this might be a Matlab HDF issue). Maybe do a if exist(filename, 'file'), delete(filename);?
Could you move the read_struct function into the sw_spec2MDHisto.m file (I don't think it's used anywhere else and I would rather it not be in the root folder as a separate m-file).
Could you add an example to the doc-text / help? It's a bit confusing at the moment. In particular, I would change the text to be something like:

% spectra: a structure calculated by sw_egrid
%     Only a single q-direction can be given in the call to spinwave
% proj: a 3x3 matrix where each column is a vector defining the orientation of the view.
%     One of the vectors must be along the q-axis used in the spinwave call 
% dproj: is a 3 vector that is the bin size in each direction
% filename: is the name of the nexus file.
%
% Example:
% q0 = [0 0 0];
% qdir = [1 0 0];
% spec = sw_egrid(spinwave(sw_model('triAF', 1), {q0 q0+qdir 100}))
% proj = [qdir(:) [0 1 0]' [0 0 1]'];
% dproj = [0.05, 0.05, 0.05];
% sw_spec2MDHisto(sp, proj, dproj, 'testmdh.nxs');
%
% Note that: 
% (1) In the call to `spinwave`, only one q-direction may be specified
%    e.g. the HKL specifier must be of the form {q0 q0+qdir nsteps}
% (2) one column in the `proj` matrix must be the q-direction used in 
%    `spinwave` (e.g. `qdir`).

Otherwise I tested with the following code:

swo = sw_model('triAF', 1);
sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 0 0] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1 0 0; 0 1 0; 0 0 1], [0.01, 0.01, 0.01], 'testmdh.nxs');

And I could read the resulting testmdh.nxs file in Mantid 6.6, where it gives a 2D spectra along the Qh direction and the MDH files says that it's integrated between -0.005 and +0.005 along the Qk and Ql directions (which bin size is given by the dproj input I guess?).

However, when I use:

sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 0 0] [0 1 0] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1 0 0; 0 1 0; 0 0 1], [0.01, 0.01, 0.01], 'testmdh.nxs');

It gives the wrong spectrum - it plots the whole thing as if it was only along Qk.

And when I use:

sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 1 1] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1 0 0; 0 1 0; 0 0 1], [0.01, 0.01, 0.01], 'testmdh.nxs');

It gives an error, as does:

sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 1 1] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1 0 0; 1 1 0; 1 0 1], [0.01, 0.01, 0.01], 'testmdh.nxs');

So - does proj have to be a set of orthogonal axes?

In future is it possible to automatically calculate proj from the spectrum structure? (Maybe if the user gives one other q-direction / defines a scattering plane).

granrothge · 2023-09-01T18:25:52Z

@mducle All Good points. I will update accordingly.

granrothge · 2023-09-01T22:38:54Z

@mducle I think I have changed things to match the requested changed. Let me know if anything else is needed.

RichardWaiteSTFC · 2023-09-05T10:16:19Z

Thanks for this @granrothge!
I just had a quick test and i think I still see the same issues as @mducle .

With the slice going along two different symmetry directions

sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 0 0] [0 1 0] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1 0 0; 0 1 0; 0 0 1], [0.01, 0.01, 0.01], 'testmdh.nxs');

it produces this workspace in mantid

i.e. it has 200 bins all along K.
I actually think with the limitations of MDHisto workspaces in mantid that we can't do better than that...do you have a requirement to support such slices? If not maybe we only support slices along a single Q direction? Or if you do need to support these slices we could print a warning to the user so they know what to expect in the .nxs file?

The following code

sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 1 1] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1 0 0; 0 1 0; 0 0 1], [0.01, 0.01, 0.01], 'testmdh3.nxs');

produces this error

Error using reshape
Number of elements must not change. Use [] as one of the size
inputs to automatically calculate the appropriate size for that
dimension.

Error in sw_spec2MDHisto (line 81)
signal = reshape(dat,Dszs);

RichardWaiteSTFC · 2023-09-05T10:20:36Z

I was wondering does the proj matrix need to have the normalised unit vectors? Or maybe I'm mis-interpreting it...?
When I specify a slice along [110] from the doc string I gather 110 needs to be in the proj matrix, but when I try

sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 1 0] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1 1 0; -1 1 0; 0, 0, 1], [0.01, 0.01, 0.01], 'testmdh3.nxs');

I get the same error as above.
If I normalise the vectors in proj like so

sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 1 0] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1/sqrt(2) 1/sqrt(2) 0; -1/sqrt(2) 1/sqrt(2) 0; 0, 0, 1], [0.01, 0.01, 0.01], 'testmdh3.nxs');

I don't get the error, but the algorithm seems to hang and the .nxs file is not closed (so can't be opened in mantid).
Am I doing something wrong? If so maybe some more info in the doc string and a more informative error would be helpful!

granrothge · 2023-09-05T17:42:26Z

@RichardWaiteSTFC @mducle OK Let me do some more rigorous testing. As to question 1. I Do not envision doing more than one direction at a time. I will at least add some documentation and may add an error catch.

granrothge · 2023-09-06T16:41:00Z

@mducle @RichardWaiteSTFC I think I have addressed the off axis case. Good catch. I also added a prototype of a test script that provides which directions should be tested. I'm not familiar with the testing setup so it will likely need to be modified to actually run.

RichardWaiteSTFC · 2023-09-07T16:16:28Z

Thanks for the changes - I think there may still be a bug with determining the x-axis of the slice from the proj matrix provided.
If I make this cut along [110]

swo = sw_model('triAF', 1);
sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 1 0] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1 1 0; -1 1 0; 0, 0, 1], [0.01, 0.01, 0.01], 'testmdh.nxs');

And plot it in mantid I get a cut along [1-10]

I think there may also be an issue with the normalisation of the basis vectors - or at least if you supply a normalised proj matrix when the argument for sw_egrid isn't normalised like so

sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 1 0] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1/sqrt(2) 1/sqrt(2) 0; -1/sqrt(2) 1/sqrt(2) 0; 0, 0, 1], [0.01, 0.01, 0.01], 'testmdh3.nxs');
sw_plotspec(sp)

you get this

But of course the Bragg peak isn't at (0.7,0.7.0) it's at (1,1,0)...

RichardWaiteSTFC

Unit test has a few bugs - but thanks for writing them!

RichardWaiteSTFC · 2023-09-07T16:19:23Z

+sw_tests/+unit_tests/unittest_spinw_spec2MDHisto.m

+        function test_1_0_0(testCase)
+            q0 = [0 0 0];
+            qdir = [1 0 0];
+            spec = sw_egrid(spinwave(sw_model('triAF', 1), {q0 q0+qdir nsteps}));


The spinw model appears to be the same used for all tests. You can add this to a TestClassSetup or TestMethodSetup like so

spinw/+sw_tests/+unit_tests/unittest_spinw_optmagk.m

Lines 16 to 29 in 49880bc

methods (TestClassSetup)

function set_seed(testCase)

testCase.orig_rng_state = rng;

rng('default');

end

end

methods (TestMethodSetup)

function setup_chain_model(testCase)

testCase.swobj = spinw();

testCase.swobj.genlattice('lat_const', [3 8 8])

testCase.swobj.addatom('r',[0; 0; 0],'S',1)

testCase.swobj.gencoupling();

end

end

In this case TestClassSetup is probably best as it isn't changed in any of the tests.

RichardWaiteSTFC · 2023-09-07T16:22:06Z

swfiles/sw_spec2MDHisto.m

+%    `spinwave` (e.g. `qdir`).
+
+
+if nargin==0


We typically add a unit test for this like so

spinw/+sw_tests/+unit_tests/unittest_spinw_addatom.m

Lines 33 to 40 in 49880bc

function test_no_input_calls_help(testCase)

% Tests that if call addmatrix with no input, it calls the help

help_function = sw_tests.utilities.mock_function('swhelp');

testCase.swobj.addatom();

testCase.assertEqual(help_function.n_calls, 1);

testCase.assertEqual(help_function.arguments, ...

{{'spinw.addatom'}});

end

RichardWaiteSTFC · 2023-09-07T16:28:22Z

+sw_tests/+unit_tests/unittest_spinw_spec2MDHisto.m

+            dproj = [norm((qdir-q0))/nsteps, 1e-6, 1e-6];
+            sw_spec2MDHisto(spec, proj, dproj, 'tmp/test112mdh.nxs');
+        end
+    end


I think there should be another end here.

RichardWaiteSTFC · 2023-09-07T16:31:54Z

+sw_tests/+unit_tests/unittest_spinw_spec2MDHisto.m

+        function test_1_0_0(testCase)
+            q0 = [0 0 0];
+            qdir = [1 0 0];
+            spec = sw_egrid(spinwave(sw_model('triAF', 1), {q0 q0+qdir nsteps}));


More importantly when I add the missing end and run the tests they fail because nsteps isn't defined!

RichardWaiteSTFC · 2023-09-07T16:37:12Z

I think there may also be an issue with the normalisation of the basis vectors - or at least if you supply a normalised proj matrix when the argument for sw_egrid isn't normalised like so
sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 1 0] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1/sqrt(2) 1/sqrt(2) 0; -1/sqrt(2) 1/sqrt(2) 0; 0, 0, 1], [0.01, 0.01, 0.01], 'testmdh3.nxs');
sw_plotspec(sp)
you get this

But of course the Bragg peak isn't at (0.7,0.7.0) it's at (1,1,0)...

Ah having looked at the doc string I think I needed to update the dproj argument - I tried this

sp = sw_egrid(sw_neutron(spinwave(swo, {[0 0 0] [1 1 0] 100})),'component','Sperp', 'Evect', linspace(0,5,200));
sw_spec2MDHisto(sp, [1/sqrt(2) 1/sqrt(2) 0; -1/sqrt(2) 1/sqrt(2) 0; 0, 0, 1], [0.01*sqrt(2), 0.01, 0.01], 'testmdh3.nxs');

but I got the same result - I would expect the end x value to be sqrt(2)

granrothge · 2023-09-07T17:05:43Z

I thought I had fixed the issue with the negative axis. I will look into it as well as the normalization issues.

granrothge · 2023-09-07T22:02:42Z

Aah here is the problem with the projection. It must be in columns not rows.

proj = [[1 1 0]' [1 -1 0]' [0 0 1]']
ans =

     1     1     0
     1    -1     0
     0     0     1

not

[1 1 0; -1 1 0; 0, 0, 1]

ans =

     1     1     0
    -1     1     0
     0     0     1

RichardWaiteSTFC · 2023-09-08T08:19:45Z

Aah here is the problem with the projection. It must be in columns not rows.

Sorry, it does say that in the doc string but I missed it! Thanks for clearing that up!

granrothge · 2023-09-08T21:34:56Z

@RichardWaiteSTFC @mducle Ok So I think this is now as ready as I can get it.

The test suite now runs and I expanded a few more cases. The files should be read back in and compared. With more understanding of the Matlab test facilities. the test code could be greatly reduced.
I now check that all of the projection vectors are orthogonal to each other (There is a test for this).
I only use the projection vector along the propagation direction to check which vector is the propagation direction. It is over written in the script so the one in the file comes from the spectra object rather than the projection matrix. Similarly the step size input for the propagation direction is over written by what is in the spectra object. I left these in for completeness. This accounts for the normalization issue of the x axis.
the CI is failing because it does not know where to write the files.

mducle

Thanks Garrett, could you change the test file to this:

classdef unittest_spinw_spec2MDHisto < sw_tests.unit_tests.unittest_super

    properties
        swModel = [];
        tmpdir = '';
        testfilename = '';
        nsteps = {100};
    end

    properties (TestParameter)
        testpars = struct(...
            'test_1_0_0', struct('q0', [0 0 0], 'qdir', [1 0 0], 'proj', [[1 0 0]' [0 1 0]' [0 0 1]'], 'nxs', 'test100mdh.nxs'), ...
            'test_1_1_0', struct('q0', [0 0 0], 'qdir', [1 1 0], 'proj', [[1 1 0]' [1 -1 0]' [0 0 1]'], 'nxs', 'test110mdh.nxs'), ...
            'test_1_1_1', struct('q0', [0 0 0], 'qdir', [1 1 1], 'proj', [[1 1 1]' [1 -1 0]' [1 1 -2]'], 'nxs', 'test111mdh.nxs'), ...
            'test_1_1_2', struct('q0', [0 0 2], 'qdir', [1 1 0], 'proj', [[1 1 0]' [1 -1 0]' [0 0 1]'], 'nxs', 'test112mdh.nxs'), ...
            'test_1_1_2_2', struct('q0', [0 0 2], 'qdir', [1 1 0], 'proj', [[1 -1 0]' [1 1 0]' [0 0 1]'], 'nxs', 'test112_2mdh.nxs'), ...
            'test_2_2_2', struct('q0', [2 2 2], 'qdir', [1 1 0], 'proj', [[1 1 0]' [1 -1 0]' [0 0 1]'], 'nxs', 'test222mdh.nxs'));
    end

    methods (TestClassSetup)
        function setup_model(testCase)
            testCase.swModel = sw_model('triAF', 1);
        end
        function setup_tempdir(testCase)
            testCase.tmpdir = tempdir;
        end
    end

    methods (TestMethodTeardown)
        function remove_tmpdir(testCase)
            delete(testCase.testfilename);
        end
    end

    methods (Test)
        function test_qdirs(testCase, testpars)
            q0 = testpars.q0;
            qdir = testpars.qdir;
            proj = testpars.proj;
            testCase.disable_warnings('spinw:spinwave:NonPosDefHamiltonian');
            spec = sw_egrid(spinwave(testCase.swModel, {q0 q0+qdir testCase.nsteps{1}}));
            dproj = [norm((qdir-q0))/testCase.nsteps{1}, 1e-6, 1e-6];
            testCase.testfilename = fullfile(testCase.tmpdir, testpars.nxs);
            sw_spec2MDHisto(spec, proj, dproj, testCase.testfilename);
        end
        function test_non_ortho(testCase)
            q0 = [0 0 2];
            qdir = [1 1 0];
            spec = sw_egrid(spinwave(testCase.swModel, {q0 q0+qdir testCase.nsteps{1}}));
            proj = [qdir(:) [1 0 0]' [0 0 1]'];
            dproj = [1e-6, norm((qdir-q0))/testCase.nsteps{1}, 1e-6];
            verifyError(testCase,@() sw_spec2MDHisto(spec, proj, dproj, 'tmp/test_blank.nxs'), "read_struct:nonorthogonal")
        end
    end
end

To make it more compact and fix the filename issue.

Otherwise it looks fine to me.

granrothge · 2023-09-11T17:49:34Z

@RichardWaiteSTFC I think it is ready to approve and go.

a function to put out MDH files

d9752ca

add needed function

a6c707a

mducle requested changes Sep 1, 2023

View reviewed changes

granrothge added 3 commits September 1, 2023 14:42

moved read_struct to inside sw_spec2MDHisto.m

df28fb6

added code to overwrite existing file

a7e84c3

added more description to the help files

720147b

granrothge requested a review from mducle September 1, 2023 22:39

granrothge added 4 commits September 5, 2023 15:59

updated struct reader

d5aed00

updated to handle off axis directions

1e4618d

slight editing to comments

a586024

add a prototype test script

32cb8b7

RichardWaiteSTFC requested changes Sep 7, 2023

View reviewed changes

added more tests

99adabc

mducle requested changes Sep 9, 2023

View reviewed changes

streamlined tests

e238e26

mducle approved these changes Sep 11, 2023

View reviewed changes

RichardWaiteSTFC approved these changes Sep 12, 2023

View reviewed changes

RichardWaiteSTFC merged commit ea5ea33 into SpinW:master Sep 12, 2023
14 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add a function to save a slice in MDH format #149

Add a function to save a slice in MDH format #149

granrothge commented Aug 25, 2023

codecov-commenter commented Aug 25, 2023 •

edited

Loading

mducle left a comment

granrothge commented Sep 1, 2023

granrothge commented Sep 1, 2023

RichardWaiteSTFC commented Sep 5, 2023

RichardWaiteSTFC commented Sep 5, 2023 •

edited

Loading

granrothge commented Sep 5, 2023

granrothge commented Sep 6, 2023

RichardWaiteSTFC commented Sep 7, 2023

RichardWaiteSTFC left a comment

RichardWaiteSTFC Sep 7, 2023

RichardWaiteSTFC Sep 7, 2023

RichardWaiteSTFC Sep 7, 2023

RichardWaiteSTFC Sep 7, 2023

RichardWaiteSTFC commented Sep 7, 2023 •

edited

Loading

granrothge commented Sep 7, 2023

granrothge commented Sep 7, 2023

RichardWaiteSTFC commented Sep 8, 2023

granrothge commented Sep 8, 2023

mducle left a comment

granrothge commented Sep 11, 2023

	methods (TestClassSetup)
	function set_seed(testCase)
	testCase.orig_rng_state = rng;
	rng('default');
	end
	end
	methods (TestMethodSetup)
	function setup_chain_model(testCase)
	testCase.swobj = spinw();
	testCase.swobj.genlattice('lat_const', [3 8 8])
	testCase.swobj.addatom('r',[0; 0; 0],'S',1)
	testCase.swobj.gencoupling();
	end
	end

	function test_no_input_calls_help(testCase)
	% Tests that if call addmatrix with no input, it calls the help
	help_function = sw_tests.utilities.mock_function('swhelp');
	testCase.swobj.addatom();
	testCase.assertEqual(help_function.n_calls, 1);
	testCase.assertEqual(help_function.arguments, ...
	{{'spinw.addatom'}});
	end

Add a function to save a slice in MDH format #149

Add a function to save a slice in MDH format #149

Conversation

granrothge commented Aug 25, 2023

codecov-commenter commented Aug 25, 2023 • edited Loading

Codecov Report

mducle left a comment

Choose a reason for hiding this comment

granrothge commented Sep 1, 2023

granrothge commented Sep 1, 2023

RichardWaiteSTFC commented Sep 5, 2023

RichardWaiteSTFC commented Sep 5, 2023 • edited Loading

granrothge commented Sep 5, 2023

granrothge commented Sep 6, 2023

RichardWaiteSTFC commented Sep 7, 2023

RichardWaiteSTFC left a comment

Choose a reason for hiding this comment

RichardWaiteSTFC Sep 7, 2023

Choose a reason for hiding this comment

RichardWaiteSTFC Sep 7, 2023

Choose a reason for hiding this comment

RichardWaiteSTFC Sep 7, 2023

Choose a reason for hiding this comment

RichardWaiteSTFC Sep 7, 2023

Choose a reason for hiding this comment

RichardWaiteSTFC commented Sep 7, 2023 • edited Loading

granrothge commented Sep 7, 2023

granrothge commented Sep 7, 2023

RichardWaiteSTFC commented Sep 8, 2023

granrothge commented Sep 8, 2023

mducle left a comment

Choose a reason for hiding this comment

granrothge commented Sep 11, 2023

codecov-commenter commented Aug 25, 2023 •

edited

Loading

RichardWaiteSTFC commented Sep 5, 2023 •

edited

Loading

RichardWaiteSTFC commented Sep 7, 2023 •

edited

Loading